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ABSTRACT 


For  event  discrimination,  operational  implementation  of  a  regional  seismic  station  requires  three  sequential 
calibration  analyses.  1)  Magnitude,  distance,  and  amplitude  corrections  (MDAC)  made  to  observed  regional 
amplitudes  are  necessary  so  that  what  remains  in  the  corrected  amplitude  is  mostly  information  about  the  seismic 
source-type.  Corrected  amplitudes  can  be  used  in  ratios  to  discriminate  between  earthquakes  and  explosions. 
Calibration  of  MDAC  can  be  accomplished  with  empirical  Bayes  estimation,  which  naturally  provides  metrics  to 
determine  when  adequate  calibration  data  have  been  acquired,  and  provides  statistical  assurance  that  the  errors 
associated  with  MDAC  calibration  are  negligible  in  future  operational  discrimination  analysis.  2)  MD AC-corrected 
amplitudes  can  then  be  used  in  ratios  to  discriminate  between  earthquakes  and  explosions.  However,  there  remain 
source  effects  such  as  those  due  to  depth,  focal  mechanism,  local  material  property  and  apparent  stress  variability 
that  cannot  easily  be  determined  and  applied  as  amplitude  corrections.  We  have  developed  a  mathematical  model  to 
capture  these  near  source  effects  as  random  (unknown)  giving  an  error  partition  of  three  sources:  model  inadequacy, 
station  noise  and  amplitude  correlation.  This  mathematical  model  is  the  basis  for  a  general  multi-station  regional 
discriminant.  Calibration  analysis  for  the  standard  error  of  the  discriminant  includes  the  calculation  of  the  variances 
of  model  inadequacy  and  station  noise,  and  amplitude  correlation.  3)  Likelihood-based  seismic  event  identification 
analysis  with  MDAC  discriminants  requires  estimated  source  population  means  and  covariance  matrices  for  the 
discriminants  from  each  of  the  possible  source  types  used  in  our  analysis  (e.g.,  deep  earthquake,  shallow  earthquake, 
and  explosion).  Anderson  et  al.  (2007)  note  that  source  population  covariance  matrices  and  the  pooled  covariance 
matrix  are  best  estimated  element-wise  to  fully  utilize  available  calibration  events.  We  propose  an  algorithm  that 
may  be  used  to  mildly  adjust  an  element-wise  covariance  matrix  to  ensure  positive  definiteness  and  non-singular 
behavior.  The  algorithm  uses  the  Frobenius  norm  as  the  calibration  metric  because  it  minimally  adjusts  the  variance 
terms  of  an  element-wise  covariance  relative  to  the  off-diagonal  covariance  terms  to  achieve  a  stable,  positive  semi- 
definite  covariance  matrix.  Importantly,  we  show  that  it  is  not  necessary  to  propagate  the  errors  of  MDAC  parameter 
estimates  to  final  event  identification  analysis. 
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OBJECTIVES 


Our  objective  is  to  develop  efficient  methods  to  calibrate  a  regional  seismic  station  for  event  discrimination.  The 
formulation  of  the  discriminant  from  a  station  requires  amplitude  corrections  based  on  MDAC,  the  correct 
calculation  of  the  discriminant  standard  error,  and  the  estimation  of  a  stable  covariance  matrix  that  describes  the 
discriminants  mathematical  relationship  with  other  discriminants. 

Corrections  to  regional  phases  to  account  for  magnitude,  distance,  and  amplitude  are  necessary  so  that  what  remains 
is  fundamentally  information  about  seismic  source  type.  MDAC  requires  calibration  data  for  a  station,  which  may 
not  be  available  for  a  time  after  it  becomes  operational.  In  order  to  incorporate  the  station  into  event  identification 
processing,  MDAC  calibration  can  be  performed  using  an  empirical  Bayes  approach.  The  approach  allows  the 
station  to  be  quickly  assimilated  into  the  network  and  also  provides  natural  metrics  for  determining  when  adequate 
calibration  data  have  been  acquired. 

Once  amplitudes  have  been  corrected  with  MDAC,  there  are  still  near  source  effects  which  have  not  been  accounted 
for  in  the  model,  such  as  effects  due  to  depth,  focal  mechanism,  local  material  properties,  or  apparent  stress 
variability.  The  uncertainty  introduced  by  these  effects  needs  to  be  accounted  for  in  the  statistical  formulation  of  a 
regional  amplitude  discriminant.  The  mathematical  model  presented  below  partitions  total  uncertainty  into  three 
sources:  amplitude  correlation,  model  inadequacy,  and  station  noise.  The  model  inadequacy  partition  represents  the 
near-source  effects  and  is  handled  statistically  as  a  random  component  in  the  model. 

The  event  classification  matrix  (ECM)  [Anderson  et  al.  (2007)]  implements  likelihood  based  discrimination,  which 
requires  an  estimate  of  population  means  and  covariance  matrices  for  each  source  type  population  of  interest.  In 
order  to  take  advantage  of  all  possible  calibration  data,  covariance  matrices  may  be  estimated  element-wise,  which 
does  not  guarantee  numerical  stability  for  ECM  calculations.  We  present  an  algorithm  that  adjusts  the  covariance 
matrix  in  a  way  to  ensure  positive  semi-definiteness  while  minimally  affecting  the  variance  terms  in  the  matrix. 

RESEARCH  ACCOMPLISHED 
MDAC  Calibration 

The  ratio  of  regional  P  and  S  wave  amplitude  measurements  at  high  frequencies  can  discriminate  between 
earthquakes  and  explosions  (e.g.,  Walter  et  al.  (1995);  Taylor  (1996);  Bottone  et  al.  (2002)).  An  issue  with  using 
these  amplitudes  in  a  practical  application  is  how  to  remove  the  effects  due  to  path,  site  and  magnitude  to  emphasize 
the  source  differences.  In  Taylor  and  Hartse  (1998),  Taylor  et  al.  (2002)  and  Walter  and  Taylor  (2002)  the  MDAC 
technique  corrects  each  regional  phase  (e.g.,  Pn,  Pg,  Sn,  Lg)  amplitude  as  a  function  of  frequency  in  an  attempt  to 
make  amplitudes  independent  of  distance,  magnitude,  and  station.  MDAC  is  a  simple  physically  based  model  that 
accounts  for  propagation  effects  such  as  geometrical  spreading  and  Q,  and  corrects  observed  amplitudes  assuming 
the  scaling  of  an  earthquake  spectral  model  developed  by  Bmne  (1970).  The  idea  of  using  an  earthquake  MDAC 
model  to  correct  amplitudes  is  that  spectra  from  an  explosion  will  exhibit  a  poor  fit  to  the  model,  which  will  be 
apparent  in  an  observed  discriminant.  Because  of  complex  explosion  source  phenomenology  it  is  not  necessarily 
obvious  which  combinations  of  regional  phases  will  best  separate  earthquake  and  explosion  populations.  The 
MDAC  technique  allows  the  formulation  of  any  combination  of  regional  phases  in  any  frequency  band,  so  that  a 
diversity  of  discriminants  can  be  explored.  The  MDAC  model  partitions  regional  seismic  spectra  into  component 
parts.  The  instrument-corrected  regional  phase  spectra  can  be  thought  of  as  a  convolution  between  the  source-type 
and  the  path.  In  the  frequency  domain,  this  can  be  mathematically  represented  as 

A(a>,A)  =  S(a>)G(A)P(co)B(co,A)  ,  (1) 

where  S  is  the  source  spectrum,  G  is  geometrical  spreading,  P  is  the  frequency-dependent  site  effect,  and  B  is  the 
anelastic  attenuation  with  function  arguments  epicentral  distance  A  and  angular  frequency  to.  Here  we  have  split  the 
path  effect  into  three  components:  1)  a  frequency  independent  geometrical  spreading  component,  2)  a 
range-independent  and  frequency-dependent  site  effect,  and  3)  an  anelastic  attenuation  component.  The  logarithm  of 
both  sides  of  Equation  1  gives 

log^4((U,A)  =  logS(®)  +  logG(A)  +  log/’(o)  +  logi?(®,A).  (2) 

To  remove  distance  and  magnitude  trends  in  the  data,  we  correct  the  observed  spectrum  log  Ao((0,A)  so  that 
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Y  =  log  Ac(a,A)  =  logd0(®,A)-logd(®,A) ,  (3) 

where  Ac(ro,A)  is  the  corrected  spectrum.  Equation  3  is  used  to  calculate  corrected  MDAC  amplitudes,  denoted  as 
Y,  that  are  then  used  to  construct  discriminants.  Specifically,  from  Equation  3,  the  corrected  amplitude  Y  is  a  log 
observed  amplitude  minus  MDAC — Y  is  the  MDAC  residual. 

Referencing  the  development  in  Taylor  and  Hartse  (1998),  Taylor  et  al.  (2002),  and  Walter  and  Taylor  (2002), 
before  MDAC  is  used  in  an  operational  setting,  calibration  earthquakes  determine  an  average  station  site  effect  (©) 
[=unitless],  an  attenuation  parameter  (y)[=unitless],  and  the  average  stress  drop  parameter  (a)  [=Passcals],  which  are 
embedded  in  source  spectrum  (S)  and  anelastic  attenuation  (B)  terms.  It  is  important  to  recognize  that  the 
fundamental  sources  of  error  for  MDAC  corrections  are  the  model  inadequacy  affecting  all  stations  and  individual 
station  noise.  Once  the  parameters  are  estimated  from  calibration  earthquakes,  the  MDAC  equation  is  treated  as  a 
known  physical  correction. 

From  equation  3,  log  Ac((o,A)  can  reasonably  be  modeled  as  Gaussian  with  zero  mean  and  variance  i2.  Here  x 
represents  model  uncertainty  and  station  noise  combined  into  one  term,  and  only  for  earthquakes.  The  likelihood  of 
n  calibration  earthquakes  is 

/[data|0,  y,  a,  r]  =  f[  fl  /[log  ACj  (coi ,  A),  M0.  I©,  y,  a,  r],  (4) 

j  i= i 

where  j  indexes  frequency,  i  indexes  event,  and  Moi  is  the  event  moment.  In  application  reasonable  estimates  and 
bounds  can  be  placed  on  ©,  y,  and  a  from  geophysical  knowledge.  The  Bayesian  formulation  of  these  estimates  are 
simply  represented  as  prior  probability  density  functions  (PDFs)  f(0),  f(y),  and  t'(a).  Here,  use  a  Uniform)/,  u)  PDF 
in  each  case  with  /  and  u  specified  from  geophysical  knowledge.  The  prior  PDF  f(x)  is  also  modeled  as  Uniform. 
More  sophisticated,  physically  based  prior  for  MDAC  parameters  could  be  developed  with  further  research.  The  full 
likelihood  is  now 


/[©,  7,  <r,  U  data]  =  /(©)/  (y)f' ( a)f  (r)/(data|©,  y,  a,  r)  (5) 

and  the  posterior  PDF  is 

/[©,  Y,  (7,  U  | data]  =  c(data )/(©)/(/)/(cr)/i (r)/(data|©,  y,  a,  r) ,  (6) 

where  c(data)  is  a  constant  that  ensures  integration  to  unity.  The  uniform  priors  simply  ensure  that  an  MDAC 
parameter’s  range  of  possible  values  agree  with  physical  basis.  Several  possible  values  for  0,  y,  q  and  x  can  be 
calculated,  each  derived  from  Equation  (6),  the  posterior  PDF.  We  recommend  the  mode  of  the  posterior  PDF  as 
MDAC  calibration  values  which  are  conceptually  values  of  0,  y,  a,  and  x  that  are  the  most  probable  given  the 
earthquake  calibration  data. 

The  data  analyzed  as  an  example  are  events  at  the  Nevada  Test  Site  (NTS)  observed  with  combinations  of  four 
seismic  stations:  Kana,  Utah,  (KNB);  Elko,  Nevada  (ELK),  Landers,  California  (LAC)  and  Columbia  College, 
California  (CMB).  Observed  amplitudes  are  RMS  measurements  converted  to  pseudo  spectra  by  application  of 
Parseval’s  theorem  [see  Appendix  B  of  Taylor  et  al.  (2002)]  with  a  6  to  8  Hertz  filter  window  (so  that  the  number  of 
frequencies  is  j  =  1).  KNB  observed  43  earthquakes;  these  events  are  used  to  illustrate  the  Bayesian  calculations  of 
MDAC  parameters  for  the  Lg  phase.  Posterior  PDFs  are  shown  in  Figure  1 .  Note  the  small  amount  of  variability 
(spread)  in  the  PDFs  for  ©  and  y.  The  slightly  higher  variability  in  the  PDF  for  a  is  because  the  spectrum  in  this 
example  was  constrained  by  high-frequency  amplitudes  and  the  event  moment.  With  lower  frequency  amplitudes 
included  in  the  analysis,  we  expect  the  variability  of  the  posterior  PDF  for  ct  to  be  similar  to  that  for  0  and  y.  This 
will  be  demonstrated  in  a  follow-on  comprehensive  analysis. 
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Figure  1.  Posterior  PDFs  for  MDAC  parameters.  The  posterior  mode  is  used  for  the  calibration  value. 


Discriminant  Errors:  MDAC  Model  Inadequacy  and  Station  Noise 

The  conceptual  representation  of  the  MDAC  amplitude  model  is 

Y  =  log(corrected  amplitude)  =  #/V«(source  -  type)  +  Event  +  Noise  ,  (7) 

where  Z?/«.s(source-type)  is  a  source-type  constant.  Event  is  a  zero  mean  random  effect  with  variance  x2  that  varies 
from  event  to  event  and  represents  model  inadequacy  from  effects  such  as  depth,  focal  mechanism,  local  material 
property  and  apparent  stress  variability,  and  Noise  represents  measurement  and  ambient  noise,  also  with  zero  mean 
and  variance  a2.  Note  that  in  this  section,  calibration  analysis  partitions  the  variance  x2  in  the  MDAC  calibration 
analysis  above  into  two  variance  components  Event  (x2)  and  Noise  (a2).  The  MDAC  approach  results  in  a  Bias  term 
for  earthquakes  that  is  near  zero,  whereas  for  explosions  the  Bias  is  non-zero  indicating  discrimination  potential. 

The  error  terms  Event  and  Noise  are  modeled  as  equal  for  both  explosions  and  earthquakes;  therefore  pooled  data 
(amplitudes  from  explosions  and  earthquakes  with  their  respective  population  means  subtracted)  are  used  in  the 
calculation  of  x2  and  g2. 

The  station-averaged  corrected  amplitude  Y  =  1 T  Y /n  has  a  standard  error  x2  +  o2/n,  where  n  is  the  number  of  station 
amplitudes  used  in  the  average.  Note  that  forming  regional  discriminants  from  station-averaged  corrected 
amplitudes  exactly  parallels  the  methodology  of  the  mb  versus  Ms  discriminant  where  both  are  station-averaged 
magnitudes.  Omitting  the  term  Event  in  Equation  (7)  implies  that  the  corrected  amplitude  at  a  station  is  Bias  plus 
station  noise.  As  demonstrated  with  the  following  argument,  this  model  formulation  is  fundamentally  inconsistent 
with  the  realities  of  seismic  observation.  The  standard  error  of  Y  with  E  removed  from  Equation  (7)  is  o2/n  (x2  =  0) 
and  decreases  as  the  number  of  stations  n  observing  an  event  increases.  This  implies  that  if  enough  stations  observe 
an  event,  this  standard  error  effectively  goes  to  zero  and  the  average  corrected  amplitude  quickly  converges  to  Bias 
implying  near-perfect  discrimination  capability.  By  not  including  the  term  Event,  effects  such  as  depth,  focal 
mechanism,  local  material  property  and  apparent  stress  variability  are  not  accounted  for  in  the  theoretical  model  of 
amplitude,  and  clearly  these  effects  cannot  be  removed  by  station  averaging.  The  model  given  by  Equation  (7) 
captures  these  local  source  effects  by  admitting  that  they  cannot  be  mathematically  (theoretically)  represented. 
Treating  local  source  effects  as  a  random  effect  Event  compensates  for  them  as  a  component  in  the  standard  error  of 
a  discriminant.  Also,  the  lower  bound  of  standard  error  x2  +  o2/n  is  non-zero  and  therefore  consistent  with  realistic 
seismic  monitoring.  Another  important  property  of  this  model  is  that  a  corrected  amplitude  for  a  single  event  is 
correlated  across  stations.  The  correlation  (x2/(x2+o2))  implies  that  large  Event  adjustment  increases  correlation 
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between  stations  because  this  random  adjustment  is  applied  to  all  stations  observing  an  event,  that  is,  the  stations 
stochastically  move  together.  Small  Event  adjustment  implies  the  correction  model  is  good  and  is  conceptually 
equivalent  to  stations  with  incoherent  noise.  Small  Event  adjustment  also  implies  x2  is  small  and  the  standard  error 
of  Y  is  reduced  further  through  station  averaging. 

The  data  used  to  illustrate  the  calculation  of  the  variance  components  x2  and  ct2  are  events  at  and  surrounding  the 
Nevada  Test  Site  (NTS).  Events  were  observed  with  combinations  of  four  seismic  stations:  Kanab,  Utah  (KNB); 
Elko,  Nevada  (ELK),  Landers,  California  (LAC)  and  Columbia  College,  California  (CMB).  MDAC  amplitudes  from 
these  stations  are  Pg  and  Lg  amplitudes  from  pseudo  spectral  measurements  with  a  6  to  8  Hertz  filter  window.  After 
applying  data  quality  metrics  (e.g.,  signal  to  noise  and  removal  of  events  within  100  kilometers  of  a  station),  the  data 
table  consisted  of  41  earthquakes  (EQ)  and  159  explosions  (EX)  for  a  total  of  200  events.  Moment  magnitudes 
ranged  from  2.6  to  6. 1 .  Amplitude  corrections  for  discrimination  remove  the  effects  of  magnitude,  source  scaling 
and  distance  so  that  what  remains  in  the  corrected  amplitude  is  fundamentally  information  about  source  type.  Figure 
2  demonstrates  the  removal  of  the  effect  of  moment  magnitude  from  the  Pg  and  Lg  amplitudes  with  MDAC.  Also, 
Figure  2  shows  the  explosions  and  earthquake  amplitudes  before  and  after  removal  of  population  means  -  these 
centered  data  are  used  to  calculation  the  variance  components  x2  and  a2,  and  additionally  the  correlation  between 
amplitudes  in  a  discriminant.  The  fit  residuals  are  given  in  Figure  3  and  the  calculated  variance  components  are 
given  Table  1. 
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Figure  2.  Scatter  plots  of  MDAC  corrected  amplitudes  Lg  and  Pg  versus  moment  magnitude  Mw  for 

earthquakes  (dark)  and  explosions  (light).  MDAC  corrects  earthquakes  to  zero  mean.  Also,  scatter 
plots  of  the  MDAC  corrected  amplitudes  Lg  versus  Pg  for  earthquakes  (dark)  and  explosions 
(light).  With  the  explosions  mean  centered,  the  pooled  data  are  used  to  calculate  the  variance 
components  t2  and  c2  for  both  populations. 


539 


2008  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


MDAC  Pg  Amplitude 


*•».  f  * 

* 

•  •  "*  •  j.  •  ■ 

• 

t  • 

i***  .•  • 

few  ' 

.i  r 
•  ; 

— 

.  .* 

• 

Al' 

-  * 

»•  •*>V.  '  *  • 

•  • 

-1  -0.5  0  0.5  1 

Event 


MDAC  Lg  Amplitude 


•• 

• 

X 

• 

• 

• 

• 

•  •!  J 
tj*  t  • 

"'W&t 
•  % 

• 

• 

•• 

*S9b£Ii' 

•  V 

•  .  • 

• 

-I  -0.5  0  0.5  1 

Event 


Figure  3.  Fitted  model  residuals  for  Event  and  Noise  using  pooled  explosions  and  earthquake  data.  Residuals 
indicate  reasonable  agreement  with  model  assumptions  of  uncorrelated  terms  Event  and  Noise. 


Table  1.  Calculated  variance  components  r2  and  a2.  From  the  pooled  data  presented  in  Figure  2,  the 
calculated  correlation  between  Pg  and  Eg  is  0.95. 


Event  (i2) 

Noise  (ct2) 

Pg 

0.23 

0.04 

Lg 

0.16 

0.02 

Structured  Covariance  Matrix  of  Discriminants 


Seismic  event  identification  using  regularized  discrimination  requires  estimated  source  population  means  and 
covariance  matrices  for  the  discriminants  from  each  of  the  possible  source  types  (e.g.,  deep  earthquake,  shallow 
earthquake,  or  explosion).  Anderson  et  al.  (2007)  notes  that  the  estimated  source  population  covariance  matrices 
denoted  Sk,  k  =  1,  2...  K,  and  the  pooled  covariance  matrix,  S0,  are  estimated  element-wise  in  order  to  take 
advantage  of  all  available  calibration  data.  Due  to  few  calibration  events  or  strongly  correlated  discriminants,  one  or 
more  Sk  and/or  S0  may  be  singular.  We  present  an  algorithm  inspired  by  Shaw  and  Geyer  (1997)  that  may  be  used  to 
adjust  elements  of  singular  covariance  matrix  estimates  prior  to  regularized  discriminant  analysis  (RDA). 


Anderson  et  al.  (2007)  uses  standardized  discriminants  Y  to  solve  the  seismic  identification  problem.  We  assume 
that  the  k  populations  are  independent  and,  for  the  kth  source  population,  Y  ~  MVNp  (pk,  Zk);  where  pk  and  Zk  are 
the  source  mean  vector  and  covariance  matrix,  respectively,  and  p  discriminants  are  used.  The  log  likelihood 
function  for  pk  and  Zk,  given  the  data  y  is 


l  Gu*.E*l y) =-^g\27&k\- ^  tr  (vX)-  ^-(y  -  Mkf  Sj1  (y  -  //*-)  0) 


where  tr(-)  denotes  the  trace  operator  and 


5  = 


X(. Pi  -y\y,  -yf 


.  Maximum  likelihood  estimators  (MLEs)  for 


pk  and  Zk  are  yk  and  Sk,  respectively.  These  estimators  are  derived  by  setting  the  partial  derivatives  of  equation  (1) 
equal  to  zero  and  then  solving  for  p  and  2.  The  framework  in  Anderson  et  al.  (2007)  can  be  carried  out  as  described 
if  Sk  is  non-singular.  If  Sk  is  singular,  however,  a  different  estimator  of  2k  is  required.  We  desire  an  estimator,  say 
S  k,  that  is  positive  semi-defnite,  as  “close”  to  Sk  as  possible,  and  that  the  variance  terms  be  minimally  different  from 
those  in  Sk  compared  to  the  covariance  terms. 
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We  choose  the  Frobenius  norm  as  the  one  to  minimize,  because  as  we  demonstrate  here,  it  is  a  norm  which  allows  us 
to  minimally  alter  the  variance  terms  of  S  compared  with  the  covariance  terms.  The  Frobenius  norm,  denoted  ||J,  is 
defined  as  the  square-root  of  the  sum  of  the  absolute  value  of  the  matrix  elements.  Our  optimization  criterion 
becomes 


minis'- (g|| 


p  p  , 

,  zz 

11*1  m'  ■ 


where  Q  is  a  p  x  p  symmetric,  invertible  matrix. 


(2) 


Because  S  is  a  symmetric  matrix,  S  =  ST,  the  eigenvalues  of  S  are  real  and  associated  eigenvectors  are  orthogonal. 
Let  the  eigen  decomposition  of  S  be  VAVT  ,  where  A  is  the  diagonal  matrix  of  ordered  eigenvalues  (L;  >  X}  for  all  i  < 
j)  and  the  columns  of  V  are  the  corresponding  eigenvectors.  Likewise,  let  Q  =  UCUT  ,  where  C  =  diag{ci  >  c2  >  ...  > 
cp}  and  U  is  the  matrix  whose  columns  are  the  corresponding  eigenvectors.  We  then  want  to  minimize 

\\VAVT -UCUT\\.  (3) 

A  linear  algebra  identity  can  be  used  to  show  that  ||fAFj  -C/C£/r  |  <  |f2||  ||A-C||.  Both  A  and  C  are  diagonal 
matrices,  where  the  diagonal  terms  are  sorted  from  highest  to  lowest  and  all  of  the  elements  are  real.  Therefore  the 
diagonal  vector  X  has  a  set  of  positive,  zero,  and  negative  values,  as  does  the  vector  c.  The  Frobenius  norm  can  be 
written  as  the  sum  of  terms  partitioned  into  positive,  zero  and  negative  terms  so  that 

I  m  n  p 

||a-c|  =  1K|a,-ci|2+  X  |a,.-c,.|2+  Z|A.-c,.f ,  (4) 

\|  /=1  i—m+l  i=n+ 1 

where  {A,b  X2, ...,  Lm}  >  0,  {  Xm+h  ...,  Ln}  =  0  and  {  Ln+1, ...,  Lp}  <  0.  By  choosing 

T,.  i  =  l,2,...,n 
0  i  =  n  + 1,  n  +  2, . . .  p 


Equation  4  becomes  ||A  -  C||  =  J  X|^;  ~ct\  and  consequently,  ||FAF7  —  L/CL/r ||  <  ||f2 

V  i=n+l 

Our  estimate  of  Sk  is  then  Q  =  VCVT  =  VDiag{  Xu  X2, ...,  Xm,  0,  0,  ...,  0}VT. 


The  following  example  shows  how  the  algorithm  works.  Figure  4  shows  the  original  eigenvalues  of  the  structured 
covariance  matrix  for  a  data  set  of  deep  earthquakes  SDeq,  which  is  singular.  Figure  5  shows  how  SDeq  is  adjusted  to 
result  in  S  deq,  a  non-singular  matrix.  The  white  areas  of  Figure  5  indicate  no  change  in  the  imposed  structure. 
Notice  that  the  diagonal  elements  are  changed  minimally,  compared  to  the  off-diagonal  elements.  This  preserves  the 
variance  structure  of  the  individual  discriminants  while  minimally  adjusting  the  covariance  terms  to  achieve 
invertibility. 
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index 


Figure  5.  Percent  adjustment  to  Sk  to  achieve  invertibility. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Three  core  conclusions  can  be  stated  from  our  analysis.  First,  MDAC  parameters  can  be  derived  with  sufficient 
precision  to  view  the  MDAC  correction  as  fixed  and  deterministic  in  event  identification  analysis  which  implies  it  is 
not  necessary  to  propagate  the  error  from  MDAC  parameter  estimation  in  event  identification  analysis.  Further, 
reasonable  physical  constraints  can  be  placed  on  MDAC  parameters  with  Bayesian  theory  to  aid  in  the  calculation  of 
MDAC  parameters.  We  believe  that  no  more  than  100  well  chosen  earthquakes  are  necessary  to  calibrate  the  MDAC 
correction  equation,  and  research  and  demonstration  analysis  for  FY09  will  verify  this  conjecture.  Second,  there  are 
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two  important  sources  of  error  when  constructing  a  discriminant  with  MDAC  correction  amplitudes  -  the  error  due 
to  MDAC  model  inadequacy  and  station  noise.  The  calibration  analysis  for  deriving  these  variance  components  is 
accomplished  by  embedding  MDAC  theory  into  a  statistical  random-effects  linear  model.  The  resulting  standard 
error  of  a  discriminant  from  this  model  correctly  reduces  station  noise  though  averaging,  but  requires  improvements 
to  physical  correction  theory  (an  improved  MDAC  model)  to  reduce  model  inadequacy  error.  Third,  the  covariance 
matrices  necessary  to  combine  a  diversity  of  discriminants  in  event  identification  analysis  can  be  calculated  with  all 
available  data  by  individually  calculating  each  element  of  a  matrix  (variance  and  covariance  elements).  However, 
these  calculations,  while  giving  a  symmetric  covariance  matrix,  do  not  guarantee  that  the  matrix  is  positive  definite. 
We  have  demonstrated  that  with  the  Frobenius  norm,  a  covariance  matrix  constructed  element-by-element  can  be 
adjusted,  with  little  impact  on  the  variance  terms,  to  achieve  positive  definiteness. 
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